{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c106181e",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# A short tutorial on USTC-Pickers (<3 min)\n",
    "* If any question, please feel free to contact me via email: `zhujun2316@mail.ustc.edu.cn`\n",
    "* Dependencies: Obspy, PyTorch, [SeisBench v0.5.0](https://github.com/seisbench/seisbench/releases/tag/v0.5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1098a5",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Load a model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4606d474",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.0\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plase specify a region or province to pick phases, e.g. Beijing. 42 pickers are available in the subfolder \"USTC-Pickers/model_list/v0.1\"TibetNE\n",
      "You are using the picker located at ../model_list/v0.1/藏东北.pt\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import seisbench\n",
    "import seisbench.models as sbm\n",
    "import os\n",
    "import glob\n",
    "import torch\n",
    "device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')\n",
    "from config import (en2cn, model_list)\n",
    "sample_rate = 50\n",
    "# 为拾取的目标区域选定合适的picker，如Sichuan（四川）、CSES（实验场）、TP（青藏高原）、China（中国）\n",
    "location = input('Plase specify a region or province to pick phases, e.g. Beijing. 42 pickers are available in the subfolder \"%s\"'%os.path.join('USTC-Pickers', 'model_list', 'v0.1'))\n",
    "if location not in en2cn:\n",
    "    exit('The region you specified is not available. Plase choose a region or province from below--------\\n%s\\n-----------------------------------------------------------------------------------------------'%(', '.join(en2cn)))\n",
    "elif location!='China':\n",
    "    model_save_path = glob.glob(os.path.join(model_list, '*'+en2cn[location]+'.pt'))[0]\n",
    "    print('You are using the picker located at %s\\n'%model_save_path)\n",
    "    # 模型初始化\n",
    "    picker = sbm.PhaseNet(sampling_rate=sample_rate)\n",
    "    # 加载模型\n",
    "    picker.load_state_dict(\n",
    "        torch.load(\n",
    "            model_save_path,\n",
    "            map_location=device).state_dict()\n",
    "    )\n",
    "else:\n",
    "    print('You are using the China picker, \"diting\", published by SeisBench')\n",
    "    picker = sbm.PhaseNet.from_pretrained('diting')\n",
    "    picker.sampling_rate=sample_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16005010",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Read waveforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "261df044",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../test_data/mseed/**.mseed\n",
      "6 Trace(s) in Stream:\n",
      "CN.BJ..BHE | 2018-07-18T15:24:11.300000Z - 2018-07-18T15:27:11.280000Z | 50.0 Hz, 9000 samples\n",
      "CN.BJ..BHN | 2018-07-18T15:24:11.300000Z - 2018-07-18T15:27:11.280000Z | 50.0 Hz, 9000 samples\n",
      "CN.BJ..BHZ | 2018-07-18T15:24:11.300000Z - 2018-07-18T15:27:11.280000Z | 50.0 Hz, 9000 samples\n",
      "CN.SC..BHE | 2010-03-14T20:34:19.500000Z - 2010-03-14T20:37:19.480000Z | 50.0 Hz, 9000 samples\n",
      "CN.SC..BHN | 2010-03-14T20:34:19.500000Z - 2010-03-14T20:37:19.480000Z | 50.0 Hz, 9000 samples\n",
      "CN.SC..BHZ | 2010-03-14T20:34:19.500000Z - 2010-03-14T20:37:19.480000Z | 50.0 Hz, 9000 samples\n"
     ]
    }
   ],
   "source": [
    "from obspy import read\n",
    "from config import (sac, mseed)\n",
    "\n",
    "# Read the 3-component waveforms as input\n",
    "#inputfile = mseed.replace('Beijing', 'Sichuan')\n",
    "inputfile = mseed.replace('Beijing', '*')\n",
    "#inputfile = sac.replace('Beijing', 'Sichuan')\n",
    "#inputfile = sac.replace('Beijing', '*')\n",
    "print(inputfile)\n",
    "stream = read(inputfile)\n",
    "print(stream)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d40ae72",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model response (optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a096abdc",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABTVElEQVR4nO2dd3xUVfbAvyeTXoDQOwk91IiIoqKiIqCiYsXK/nQXse26a6Ooix11d107oqLYO8KiiAXBhlTpRWlKpIP01Jn7++POZObNTEJCJiQZzvfzmc+8d8t7Z968986959x7rhhjUBRFURQfMVUtgKIoilK9UMWgKIqiOFDFoCiKojhQxaAoiqI4UMWgKIqiOIitagEqSv369U1GRkZVi6EoilKjWLBgwQ5jTINweTVeMWRkZDB//vyqFkNRFKVGISK/lpSnpiRFURTFgSoGRVEUxYEqBiUyzH4WJg6qaikURYkANd7HoFQTpo+qagkUpdwUFhaSk5NDXl5eVYtSaSQmJtK8eXPi4uLKXCciikFEJgDnAtuMMV28aXWBd4EMYANwqTHmD2/eSOA6wA381Rgz3Zt+LPAqkAR8CvzNaDAnRVEqiZycHNLS0sjIyEBEqlqciGOMYefOneTk5JCZmVnmepEyJb0KDAhKGwF8ZYxpB3zl3UdEOgFDgM7eOs+JiMtb53lgGNDO+wk+pqIoSsTIy8ujXr16UakUAESEevXqlbtHFBHFYIz5BtgVlHw+MNG7PRG4ICD9HWNMvjFmPbAG6CUiTYBaxpjZ3l7CawF1FEVRKoVoVQo+Duf3VabzuZExZjOA97uhN70ZsDGgXI43rZl3Ozg9BBEZJiLzRWT+9u3bIy64cvj8+/PVVS2CoigVpCpGJYVTX6aU9NBEY8YbY3oaY3o2aBB24p5SRTw9Y01Vi6AoNQqXy0V2djbdu3enR48e/PDDDwBs2LCBLl26OMqOGTOGf/3rXwD86U9/IjMzk+zsbLKzsznxxBMjJlNljkraKiJNjDGbvWaibd70HKBFQLnmwCZvevMw6YqiKFFLUlISixYtAmD69OmMHDmSWbNmlanu448/zsUXXxxxmSqzxzAFGOrdHgpMDkgfIiIJIpKJdTLP9Zqb9onICWKNYtcE1FEURYl69u7dS3p6elWLEbHhqm8DpwH1RSQH+CcwFnhPRK4DfgMuATDGLBeR94AVQBFwkzHG7T3UDfiHq07zfhRFUSqd+/63nBWb9kb0mJ2a1uKfgzqXWiY3N5fs7Gzy8vLYvHkzM2bMKM5bu3Yt2dnZxftbtmzh9ttvL96/4447ePDBBwHo3Lkzb775ZkTkjohiMMZcXkLWGSWUfwh4KEz6fKBLaA1FUZToJNCUNHv2bK655hqWLVsGQJs2bYrzwPoYAqksU5LOfFYURYFDtuyPBL1792bHjh1U9WhLVQxKxXj6WGiv8xAVJRKsWrUKt9tNvXr1OHjwYJXJoYpBqRg718DsZ6paCkWpsfh8DGBDWEycOBGXy1V6JS+BPgaAuXPnEh8fX2GZVDEoEeXHhJugqB/EVvzmVJSjAbfbHTY9IyOj2NfgI9DH8Oqrr1aaTBp2W4kojeUPOLijqsVQFKUCqGJQIk9RflVLoChKBVDFoESep7Lhp8iMp1YU5cijikGpHFZNhfz9kPtHVUuiKEo5UcWgVA7GwH+7wqMZVS2JoijlRBWDUnnkepfoeOdKcBdVrSyKopQZVQxK5fBzQJirVVNh5y9VJ4uiVGMqEna7WbNm5OfbwR47duwgIyMjIjKpYlCODLp0t6KExRcrafHixTzyyCOMHDmyzHVdLhcTJkyIuEyqGI4WqtqUs3kR7NXlNRSlNMobdvvWW2/liSeeoKgoss+3znw+Glj8LkwaBn/9Ceq2rvjxjIFZj0H3IWWv8/EN4IqHe3QpVqWaMm0EbFka2WM27goDx5ZapCJht1u2bMnJJ5/M66+/zqBBgyImtiqGo4FP/mG/t66IjGL4Yz3MfJiVM14nqzx9TndBxc+tKFFGRcJuA4waNYrzzjuPc845J2IyqWI4GijYH9njeTwAJKIveiWKOETL/khwOGG327ZtS3Z2Nu+9917E5Kh2PgYRGSAiq0VkjYiMqGp5ooswDuADO+GT26GoHC/5nWsAyIzZGiG5FEUBZ9jt8jB69Oji0UqRoFr1GETEBTwL9ANygHkiMsUYs6JqJYtiPh8Ni9+GFr2g26VlqzP11koVSVGOJioSdttH586d6dGjBwsXLoyITNVKMQC9gDXGmHUAIvIOcD52fWglkhQcgLhkcBd6E6TsdXXoqaJEjEiF3f7oo48iJlN1MyU1AzYG7Od40xyIyDARmS8i86t6CbwahbG+AXb/Bg83hbkv+tOkHIph/5aKy1JUAD88A56Ah+K5E+GnNyp+bEVRKkR1Uwzh3k4hzVNjzHhjTE9jTM8GDRocAbGihIWv2e9d6+z3zEcovrw586EwF8bUhsXvVL4sDzawZqwXToH138LmxbBtOUy+qfLPrShKqVQ3xZADtAjYbw7orKhIsWu9/faZj3J3+XsMc573t9YnXe+sl78PVv6vcmTaugwmnmsVhI9pd1lTl6IoVUJ1UwzzgHYikiki8cAQYEoVyxR9+JQBwIrJ/u15L4eWXfk/eKQ5vHsVrJ0Rml9eNs6DDd+XXmbOOPjhaf9+7h86a1pRjiDVyvlsjCkSkZuB6YALmGCMWV7FYkUPh1pyc/vK0LR3r/Jv//g8tDm9YjK8fGbZygVOhvt3RyjKgzF7KnZuRVHKRHXrMWCM+dQY094Y08YY81BVyxNV5O2BF049/BnIv3wOeXsjK1NJ7A4Yg1CUZ7/z9zmd1YpSHjYtsvN2lENS7RRDpbB7o94QPjYvcpppSmPv5tC02c9GVJwSWfoerAiyIj7SHF45+8icX4k+xp8KL/atailCOJyw2263m+zsbMenfv36XHbZZRGRqVqZkiKOuwgeCJhBGGyKMKZ8wzSjhY1zDl1m76bwDue8I2jOee/q0P9s4492BFXznlbZP94aThsJp42AwjyIS6xcmZZ9BOmtoNmxlXuekjDG+lyS64bm5e8DBBJSj7hYNYbdv1a1BCEExkqaPn06I0eOZNasWaXWcblcjhhKmzdvplevXtxzzz0RkSm6ewzvXO7cD5yY9crZcF8dOzyzvBQV2KGdwRzcBcs/tmsdVzabFsGit6Aov3KO/58s2L46NH3O85VzvpJ49dzQtJfOgD058O8Odn/mIzC2FTzUCLb/XPKx9m21/3d5//Nff4DV3oWHPvg/ePF0mDm25DAin99tz3G4SnTjXBvlszAvNO+Hp+GxTOegAR+PNIdHQqb9HD4rp8L//gYFB8tf949f7TVYNzNy8lSEg7uq9vzGwKaf7HUphb1/7CS9dq1yHtowdOhQ7rjjjpAexuFS83sMm36yN+BZD8GJN9u0j2+ERW+Glp39rL/MrwEjY8adDMO/K/s5XzjFOmrv3QUx3qnruX/YB9bHvX9ATIDePbjL2co7uMvWebqHDUd9168Qn1y28+futt1isOGs79kBrriyy19W5ocZpXSk2fBt+PQnOjv383bb7wWvQt9R/lazMbYBAHDt9PDHWjkV1n4F5z5hy797FdRrA/3ut/mvDLTfIwL8HjMfsZ+7t0Fsgj/dGL+pbmxLuORVaHoMpDaCuCR/uQM77f2RFCb2/sv9/Nt3rIWU+v79L7wtwveucfam/p0V/rcdLrm74d0r7faCV53nMgb2bbb3baBsgXzzmP1+7Xy4Zye4DvNVs2UZpDa0n8Nl4esw5eZDFnt07qOs2rWqhFxjGwKu+JKtDL5glfH+HlvHuh25q9ddsMd77+Tusvdmst+S4Q+7ncvmTZuY8d644kZsSNjtTTncPvxq2L/V3lPAE088QWxsLLfcckvJP854/A2i0Vuc92IYoqfH8Plo2LnWbodTCr4yEDrscstSGFPHuZjNjjX+1uXMR/3pv83xj965v65tFRoTuuj9wlf92xvnWqXx3X9hwkCrFJ7uYT9gncEPN/G3srYsK70n8FZQTKMH6tsJYpsWhZY92py1Pz5rW83/u9Xu+5QCwIT+/u0xtWH+BPuwv3ul3X7rMlj1iV2K9PsnYfpop7N9bOAUGy8fDbPfyz6y99U3jzvz3/8TPNkdHmpcHJWWZR9ZE9ijGaHO/CXvO/cfb1Pyby0OZwLsK8Nw3sXvwCMt7f11KB5tVXLefXVsj/LxNvDR9aH5819xzmB/oJSAcLvWwapPS84fdxL8q134HjrY3zJ9dMlhWgrzQpXCrMdLKHsQPEX2EzikG+y8Gk8hFB4IL0tgBOOC/aH1Dwb4OHf/5sjymZJWzXiHz954mmv+di9m82LAH3Z70aJFLPppIcOvvshW8g7fXrx4Mf/973955ZVXkPx99n20b4ttMAdekz05/u3PDr1CnJgaHvemZ1OXmT8swKY6Zk/ppoIxe+A/nWDv7+HzT/qbbSkGH+OGH6BhJ+eLxseAsfBZmECwvlbW4ZirAE6/Gxp3h7cu8af1uMY/gzkc138LTbrZ7T9+hSe7Hd65y8lBEY7PaMHYbTs458BhmB4qg26XwZJ3K3aMy970t5wryoUv2kCFwfeD7z4pKrAzwoNJz4Cb59uXQeD/Wb893DwPVn8Gb4dxOv5jpe2RPJkdGsZkzB7IWQBJdeDT2+191Xmw7cUGN3IA/u8zaNXb2QPz8ZcZtvWa0hAObAvtzYHtWXncsGO17UF9cB20PcP2eAGu/BCaZsPKKdDuLKjdHF6/0PbkAJLqwom32GuQNQgadYZfvoCPh/vPccMPNh3sNanbGp49LlQW3+8HVs79mqzmdcKXadzNWgQ8Htiy2JnXsJOzp7jpp9D6TY+x34W5sH1V+DwgNTWV/fv22YEhQKPuZ7L0q/c4mNyccy+4yB8vadNPjPn3OFJTkrl9+DXkJjTkuFP7M3bsWM7td1rouurxqVC/HXjcrJzzJVnTAxqUY/YgIguMMT3D/fSab0oKJlAz+rhzvd/Ms+G7kpUC2Jbi90+Gpj9/ojVXhSOcUgD7Yi6pq10WZjwYmlaaUgB4oY+96TfOK/ucgQiwJdaa1MbVqV19FENFlQJETikAfPQX+D1M9MsDO0rvGfyxAbYu95sPfezw+lPCKQWwrfqSCFZOa2dYxRBOKQC8MgAkBjqGWQzmxTLMbRnXx7aUiwJa28s+8G+/eZGz/C0L/UoBrAnmq/vs9rwXw5/j+ROtApr695KtBoG4C0sfur1lCTTpDvvDjM7btgKaZFuzUkmTL/P32d8c7hy+gS9bltnehVcprFqzHrfbQ7302hzctMZf/mDoqMrbb/s7p556KueePdDKGkzBfnuew/B11XzF0LATENA1+3KMf7vLRfZmD7Ttz3vJv93mdGjXHz67q2zn8pmiSuOuX/3d8MpurV/yqjVVBHO4PZQIcBSO8Sof4Zz34ZTC35Y4759gpeDjqWPCpx8OwWYsgL6j4Wtvg8h4yhcapf0A+Pkzu70jzECG0vCZWcvLg2X0Rexabydspp9RernNi0vO27MR6rS09v5w7FwTPh2s0khIA08huXn5ZPezy+QaY5j43/tCw24HmZ82bdnOcxPfp2PHjmRndy82K3bu0IY3nwlowObucta97ecyjcSs+aaknj3N/Onvhr+RAh1mL58VOkzzzvVWaeTvL3k0R9dL7Zj6QC551Sqc8ac5u5C3/Qxpjcr3Yr7hBzvPoqRWX0mc8U/o8w87MsnXHa9C1sTFMbh5EwCWrv/tEKWVQzJmj7XVh1v7ovVp4Uf7jN5ifRmRYtgsa94Jdz/fOMf2MqaXYK8+nGehMsm+ChYF+D16DIWFE1nZ/z2yWpVRmTTIgtj4kpVFSgNrutvxC2EXxWrYGXZv8McBS2kAB4KiQ7sSIK2xf1htnVaQUAu2BqxF3TALtnn9nI27he8tBLHy121k/fEVnP1YcVpppqTocD7XC9Pi+nOQg/lPYRxcvp5EQqp9EP8RZAccvRUuCtNt7TzYfg+b6U8bNtM+CAD/Ny20zm0/23PcEmRKaNQZOgyA676Af+62duFRYbquV38Mf/oE/rrI5vfxruOcfQX0vTu0vI/sCJpCSqFKmxfhrldlcdM8uOhlOyKtNO5cDxc8bxeDv2tD+DLt+odP9zVoul8ePv+KMC17sCNN7lwfeo8dDmP2WKUAEJcSmt+wI/S+MXzdU+/yPwuR+G/uXO/fvjWM4/zm+ZBZQo8KrHnpgmfhsjesQgBYONGfn5huX8ZNsq3pqCTiEq05rWGn8PlpjSE+Jfwx6rW1SqVWQAM0WCmAva6BI9V2/+r0T6RnWOXhI1ApNOpqfRdNj3Gex8fAR0PTSiA6FEM4mgdNQCrLcLlaTezQT7DDyXyTpc4PmO177efOOmP2eB+igC59qxNDj+17UOq1geHfW2ffvX/481v0sl28Wk3tsNU711tnnu8cbfpCxslQNzN0WGuf22DEb87jgf0tFzx36N8dYQ4cyUmDN88/9DDfMXsgo0/Zj5naKPRa+qjfDrpebJ2SDQOcrNd/499u2Mk2OrKvsMOgww1JvXM9XPkepAcMcc7o4+zlxiXC0KnOejGx9gVzTdA8hosn2O/kuvYeu947zDextlVMd22w90Pg7zrxFhjyVmgjamSQn2540JDh0QGO7OFBARGHzbTDhX0E/zenjYSe11qFcdM8OO8Zf96gp2zjKJA/z7C/yfec1Wnh3/Z96reDoSXE2kxt5HcSZw2CE8Ios/RWkNbEPn8SY5/lxqWYgQOdzj6S69r/BuxxGnUBxCqEpsdYsxFYxRFM3YCGrcTY+nEB180TMPosKb1kU1DgOy4lzCCGcjyXNd/H4OPubWW3LwK0OD58uisudLZt9ytg6j+gcRdoaevlFbrpeI+1n84ddQYNawXNuL11mR2K1/w474zUABp3gfMOEZYiuS7c8UvpZXzExNgXANgHa+UU6DjIOY+ikpmb5H9YTshoUSnmpO3nvcGUj97kutiAHln9dvb70tfs2P6S+NNUh1lj418302LPfJg4CO5YByn14NfZePZt5duiLNrtzaepr3BGH+v8zd1FgdvQ/u5Pua1fe2758xcw61HbCq3XBv78lZ18d1GY+R8Bo+Xcp4zA5eut/m2RHfVS0n+V2ceO0vnF2yC53LtWRsYpcMxV/mGhnS5w1mvSreSgg+HSz3rI70PzvcR81GsDsd77++4ge3rjLvZ3r5gMZz0Q/nyn3GGH8fpMUz4atLef7pdD/l5/D/7GHyFnnh1VFu4lXBLHXA0/ve5MC7YCpGcEVTIYwvjGYlxQr51/pE9aE2d+3Tawa61/v07Q8F5XnPO3BhKf4gwrn1jL9iwDzfr12xc7pP0yBcxVatK9dP9HgBKw7oLyNdaiw8cwf77dCbRnBtz8BUUetu3Lo/msO/x2xru325ZXGIrcHgrdhqT4ktddfWbGL/zrc/8s2w1jw4zW8PLbzoM0rJVAYpz/eBkjPjlkvYhRyXber5KTuLWRs4VSGYphzY2/c+5/Pmda/AgyY7wvqMCX3C9fcP6EFUxOuNdZ0Vum48iPuSDmW5Z6MlluMln38NnExDgfmGe/XsPj062jdEPiFf76xoAxZIzyK6Xg/84YwxtzfuPrVdsYf/WxxLqCXvYeN71Gvc1x3Trx7BXhnav5RW5iY2JwBcpVmFvsO8jIewuAqbecTJdmte3cHePxK8gSKHJ7MEBcsExlxV0EmIhNpNyXV0haYiVMygQ7yuvH523PJSbMMxzwPKw/7VnSugygXr16SLgWde5uq5zCTQjL328VZinWiK1789iXV0jbhgHK1uP2m4DSW0NSCc/n3t9h/zb/fpPutkfhI9C/WaelY9Kcj6Ktq9i9Zy/7XOlktnHeI0fPcNULxtlxzT38LUffCxhgxegHSfYphth4PB7DG3N+5aS29WnTwD8Xou+/Z7Jxl3MSS/BLIFApAOw6UEDdFL+imbt+F52a1qLLP+1s2/6dG/HC1aH/QcaIT9gw9hy278snLTHWoTw8HsOWvXnUSorjYH4RH/30O5f2bOE4T0XZadKoJ/sOXbAUZiaHPjRubNz0SHLmf2YBCfQteILOsoE9JDPwkxVc2KM5WU1q4WlzJkuMP4zEdlObBqNs1PZdBwrIM3G84/YPrWw96tPi/9UYQ+ZIpx/qhLyn+XGUdxayyCH9KLe/v4QPF1ozzFtzf+Oa3hmOfCMxbCOdT5Zs5ua+e2lVL5nkeP8jaIyhw922Fzrl5pNo3SCV1IRY+1I67xlmFGTBx1Yhnvu0nalfPzWe+Xf3ozQ27jpIn8e+BsI3RGb9vJ0W6Umc/m8bn+evp7flH2d1cBZyxbJpdy7pyTHMWLWNfp0aER97eErmwakreOm79fzpxAzGnBdmzgNw+r9nsm77AX66px/p5b3fU+rDGWWLGdT8+PPJ2bKDylgiOLfAzc4Ddqjqb2uhSZ3A58TbA9uzidLWIisoisdVdABXfBLscY7sKiiMJf6At+6eRGCbI98Y+H33AXbuPkBB7XpkljIiOpgarxh27A+YIZx9OWRfzqmPf82vAQrBx98mreXFG2ZDcl1+353LSWP9ttUPhvfmmJbpzNuwK0QpgG1xhbQAA5i5ehsX9mgOwAuz1vLItFUgRbiSf8V9sA3Tl2+l0O1h7LRVvPzdekfdQOW1/pGzERGW/b6n+OEPZOy0VYwY2JE/n5zJtn353PHBYm7u244uzWqFtMC++2UHV708hzsb16Z1YSEDwswvWOhpTz/XAuZ72tMzJnycoXfSUnnMPZh6f3Tli4Q7Q/LnJYYGrvssJfmw5jNk5r1BGrksSfxLcdoGTyOmeHo7yi03GQC8+O163pm7kaX39WfF5r2YgC7zufkPsffB71n5wAB+3lq68tseeB952UI91hfUwucFGDXJ6fj0KXWAz5ZtLlYKAPdOXk6/To1oUjuJbfvyaJiWyHUT5xfnD3zS2u3XPDSw+L56bbY/js7w1xewaU8eb/3leE5sU58tbS7h2kcCxvV72bG/gBMf+YofRp7B1CWb2JdXxEU9mjte2j6lAPDBghzOz27q6DkMnTDXccynZqyhbaM0zutebEyjoMjDiWOdvoiRl+fw66Z03ppl7eaPXdyNS3uGzg4PfNYW3HMKr699HHH159UfNnBT37Y0SHOaix6cuoJ1262p5ZgHvii+xh7joftr1rH7fx1uo0vaAM7qHH4UVs8Hv+SE1nV5xtsz27Injx3782nXKJWEEb/ZcCVAXFIamZlpIfU37c7lrg+X8OyVPah1GD2b/CJ3sZL3EayUi9we2o6eFjYP4NHPVvH8zA0l5meM+IRa7CeXRFY/3DGk9/vSt+t48BPvMr6s4pqTyq4ZarwpKaFJO/Pm1K+5+Fj7Uh7w329YtaXkl8DM209j+BsLSi0TjgeucLM1/2duP+52/jhQwDEPfBFSxvfn+V70iU3fIa72IgD2rRxb5nOtemBAsf+iPJzdtTFz1//Bjv35zBt9Jsc99CUAaVl2Al6geceDtTpOKBrIdbHTeKDwSv4W+xG1xKkU98YIJ7WyD3vBrpMYtqOQ2+I+cJTpmtkyRJYL9+3nvh125M63SYkkGsNxeYcO+OczlRSbcYBueePZS9kjhnaSDfxmGrIf68BbcX9/Ot0bPk7SZT1b0L9LI5bk7OG/Xzp9OnHpswEPj591E6e2bxD2P/e1rjPCNEQAJHY38XW/46HT7uT291aELdOlWS2m3tIn6BhuJPYApqgWS8acRbcxQYMepABx5WGKbMC1Oslx7D5onZRJcS5WPjAAgCe++Jknv3L+rqa1E/lhpB2/X5Lc4HwZnf3kt6zY7A/fERO/lZQ2TwCwb+Uj+GzY953XmaEnZhSXC+6Jxdf/koQGX1K0L4vcnKEh5wknky9/3e51nD/5/OL0fSsf4fO/n0r7Rs4X+5Nf/sITX9pGzkU9mvPbrgPM2+B3ut/Utw2Z6fEM6NyQ5KSU4hfq9n35HPfQl9RLiS9u6Qee/6FPVvDit+t55opjOLebVZp5hW5iREJ6T7PX7uTyF390pL0z7AROaG3NPdv25dHrIb+ivyC7Kf8d4h/A4vEYWo/yX7cpN59Et4AZ2jv259PzwS+L96f9rQ9ZTZzB94Kv46QbT+SYlv6BEJU2XFVELhGR5SLiEZGeQXkjRWSNiKwWkf4B6ceKyFJv3lPiNeyJSIKIvOtNnyMiGWWV4/b3F5Mx4hO+X7PjkC/80/41s9xKAdcBHvtpNBNXTKTtQ2MdL4jHL3aOXjjvGX8r36cUykvZlYKHuDo/EpNghwR+unRLcQ/KpxQCWeuxDrSDInTPbMn4OrX40G1H63zpOZYLCvzOw/PyH2BIwSieTq9TnBZf93tek16OY26MDd/p/CjN/yK/sXFDrm3SiDdrpVIYtrSlAIhL/x6rtvwUlbNju8JkFCsFgOFvlDx88935G7n21fkhSiEtawSJjSeT2Ph/3D7tlbBKAeDp72Yz+n/hA/2lZY0gtd1Y4ut9x+hZjzjyxLWPtKwRxNX9hmW/7yXnj8DelZu0rNGktnsYJJ9nvw6dKJXW8V5S2z2MxFr/iU8pAOQWuun9yFdc++q8EKUgsbvZ1/RWznnvCp4M+s3BzFnnn20bqBSAYqUAEF/f/4L755TljhfSGz86o4kmNLD3ZWyaf7XAjBGf8N78jeQWuLnzg59IbPI+MYkBjRiPYcf+fK6aOtxxLInfwVlP2NFg2/bmcecHi8krdPNizsWkZY0gJn4rHy7McSgFgGe/XsvtH62kywOzaD3qU9Ztt3GOfM9MoFIA+GKFNd+9+K3t6d/81k/8vHUfl74wm473fMYpAT0ygJw/DoYoBYAh422aMaZYKcQkbCYtawSfbfsXbo+/kf7JUucw3wvGT+K6tybz1hx7Xfo+PtORP/DJb2k/eho3vVXyvb43r6jEvGAqOmxlGXAh8E1gooh0wq7X3BkYADwnIj6T8/PAMKCd9zPAm34d8Icxpi3wBFD2Qbdernwp/DoDN/Uth3EtiAt7NCO+rv9ln9TcOdX+gmP844WNMSzJKWEkiJT9TykrKe0eJrHJx6S0fpLgl6kf/3m3e/+Bvd4RMO+npbLCZJCR9xa/p25je8enWBpv7blLTBsWNVnMO7WcrTFPu2cc+2e3aEpZGVuvLg/Ud64j8H7RKfTLt5Nu+rVsRmLj/5GWNYot3pmfg/Pv46DXHitxu+wDn1C+sfHfrPnVcf3FtZ/4Bp+DhA+HEJO40bGf1OxdAq+jj9hai0lt82+m7AodAikuZ+j1+LqzHfup7e3s1MRGtlV48qP+l0tSC/8Y+4SGn/PCrHWOuoG/JbWdU+H4ZpRs3pPHjFXbCCa1ne25/pa7lCdmLA3JL5a3wXSu/eoiuo2ZzrwNpc/ZiK8buob3t79Ym/09kwNX5g26R2P8/qA7P1hC1r2fMWXzo8TVWUBK5nNYTxXc9eESej74JX/scQ71FLH5GSM+odfDX/He/Byy7nsXEXueQOVlKxRBzEFi4rfiSvmZ+AbTAcPp/57F1CUl2/n/8tp8VmxyKsaznviGuevtddmyN48f1uxgT24hj322yvFfPn5xF1I7jiSh8cfFaRO+31C8bZ9diKu9hDaj/EOTb3nb71h2Jf9CSpv/MLfwbkZNnseijbvZlx96Pxa4PXyyZDOTF/3uUDJJXr/l0Amzuff7e9m4b2NI3WAqpBiMMSuNMeHmup8PvGOMyTfGrAfWAL1EpAlQyxgz21gb1mvABQF1fE/EB8AZEnaYQCiu5DWEe3ATGn9EWtYIXttyCQmNPwrJl7gdpLS1rbJg1j18NhvGnsOgbk1JqP91SD5Ys1SgrdbhvIyxLcCsujZeTWxqeDOCv3weKW0fIS1rBK6gsjHx2zgmI/Svion1v3ySmoeLoVREWpZ/8tsvic7L6b913CQ1tyacK5r5bbZxtcK/OLpmtiy15e9jQ2xsiLqalJbK4GaN+c7TkTWeptxRNJxfTHN65T3LroAwAP1aNuP5OrX4yfhGUnhIbWsViO9hApDYP0hs+haxqeGXBo9JzCGtw/2kdfRfh+RW40ioP4O0jveGrZOSGbpKXWxA72/e6DOJid9GUrO3i9Ncyf6hi49d3I3U9mHiXHlfdMvvDwrDEOM038Wm+n09sbWcQxLfH96buy4J3/hIyxpBWtZIkjOf9D4TpZuJY1NCewwL7+kHFJFQ/2ti4ndS2OgZLhnnVGrJCc5/X1yhPrmrX57Lxl1OH1NMojNGWWLAyxIgru63xNXy/48SZ1v67y/IATzEpq6hcG9n8rfbGGCxaaGT3WISgxoNAc92SuZ/SetwPyltniC55QQS6n9NWpaduX3zW2GC4AVw9lPfElfnx6D/w1Os3K54aQ7d7/uc52auddRbkvcyIob49B+R2N3kFbp5YKrv+Xb+P/ENPmfUpKVcE+TzSW7lH/6c2PRdLnjWr4j/c2noZLq/vbOIr1b6hxV/dOOJEJNHWtZoJq2ZxNkfHXoVxMoa6N4MCFRLOd60Zt7t4HRHHWNMEbAHCBuvV0SGich8EZkfk7CJ5FYvkZZ1d7Et/axOjdgw9hzi0/0XOD59LnF1Ai+4m9S2/yImbi9pHf/pfahGEFfnR16/rlex3bFvR//cCHeutbW7vA9URv3AFkwREutvVdw22N4wp7U4DYDEgJcIwPBT2zD2wq6sf+Rsrj1/AWkdxhATZx/45BavkdrhXtY9fDZL7u9DSpv/sCbpTsDw7BU9WHhPP1ypKx3Hi01bRVrWCCRuJ20bWjNOoFIAeKxxLNlplzI92fYCtsXGElt7PmlZzhhQ3RJuCr3oQfTIbMmpLUtfFGZQi6a8nxbqG1gTH88NbQ4yuE0scXXmkNDwE3alhQY/fC69Dr5WZlyd+Y68tKwRxCTmkNruUeJqLyGpxev2/5ci/n2J/2FJyQzs4dhjxSTsKE4JNINY/A2Mfavv4+Bv19k6sdb82LhWImmJsaS0+Y+jVnIrO0P+q9tOdThgD/46DOOxprCYxC18c0df7vjmH466gS+4b+46yZEXqPwBjsuoy3NLg8NGewhsjbsSN5Pc6iUSm/hDuax7+Gym3eqcu2PNdpaP/t6EL+/qQN2UeC47y/9yjk3egLj2kdB4EkghG8aew+l9Qk1nL17jsyQbfC+8PkEmlrotrJ/kpX6vWKkL/EOcXclrSGzktIm7kvyvkKSW9vqKuHEfaOuVP9Rck9DAGXXAr7SM438PJCZhi1dmZzOmVqLfhBmT+BuJTT52NAbSskaR1mEMcXXmEpOwieTMJ4mtvQDfPXR5rxZ8vHZScfnEpu/y2uwNxfuxac6GV0L9mbw15ze++dk/Qmr+vUH/mc8EJ4XE1/+Ss7rUAez1s71p+1tum9+v+J2W1aTWoRumwdfkUAVE5EsRWRbmc35p1cKklTTLwqc2S8tzJhoz3hjT0xjTE3H+mbG1FjH+mp58vObjkHqJTWyv4dO/9iGxWfjIm4lNPiY93a9tN+61N2dSYXcK/rD29eSWL3PPuf7hfGPO60ha1t2ktnsYV4pt7RWJtc/+X5f/sz9ODPdf0IH/XNqdv57eljv7d2BIr5bkFuXy/s+hIQ4kpoC5W+fwz+//WZyWljWSc7o1oW5KPMlec8NJTZyRLVPbPs6X/ziVq04KvwqUu/6P/Ku+Py+p6QchZUzT0EBpeVsHhqTtCg70FYYHg0xHwSQ2mUR8vW9DTHQ+0rJGWXt/k9Aen/Ol7y3f8W7i6/zE1FtOJiPD+TC4kjY4lDdAQgOn7+Cp6wJGyHgScB9obcs1tM7rxy7uxth5fl9M3hb/CnMdG6fRpkEqeUW2UXBu5mDW3nsLB9ZZRRCbtoSW9ZL59nf7Yn2h3wsAxNfzv2gv/MQfEbdFmnXqS+xeklqOp3OneRws9LfCj21kZ/d/9o8eDDk1dNXAuDq2FTzxWtvQ+XKzs3ESk2BfQOef/DtDPxvK4CmDWb1rNZ9udJZLbf8Q8elzOL2X9XXM+t2+4HNzrqBwn+0R9+1Yj7UPDyAta6T3M8Lbi7OP8PpHzqYozprEejY+hrqJdRl8nH/4ZnIrf3DLgxv/VPy7fcSmWPt+3pYLcedmWPm9SjM2bZl9KSb9iivRPrtD2tp1Mi4+Yw0L7+kXYtoLJKX1f70yj+Lsro1BCknLGoHJvJ14r6KJrzfLL0utRcTW8vcwEpt8RErrp3Albiap6fvFjdSp+5zhaGJT1vPwp/4Gna+X/vDJDweUclo+bpt1WxiJDWkd7yGhwZf0frs3Yy/sUnz9Ulr/N8QUOvLbkcX3b8HOkziw7q8lXgsfh1QMxpgzjTFdwnwml1ItBwgct9YcO1g3x7sdnO6oIyKxQG2g3OvxJTV7h7yiPO753o5j7tuiL0uH+jXzlNta0KlpLeJq2Qkm/zntPyHHuPwTG6Nm3e51nD3JdrtOadeYyUNvLS7z37WD+eF3u2i3p5b/pkluOYHhp7ZmwvKXiZVYkmL9N//jqwdzTJtC/nFWh+IeyfFvOVsE/zzef5P85fO/8OVvTify8p1Ok8nTZ/6LpUOX0j39lOK0fHc+k71275OaOVugZaVR53879gt3lRKLppox6rtRXP7lyexMcprXkjPGW4duEGlZI7ymF1izxyqTe3vfy4ax59CsjrPHs2bvMj785cPi/cI/Ti7evn2w9Vmc/r5V1geLbA/wjaHWjZZQfxb7AxZ0ObGpDZ3iSrC+gMz6SeS7renjtk4vMKi1VTqp7R4mNmUdv5kPHffL2Zn23lzyx/d8ss3+riEdhpC/w/9fndE5kT5tbej38UvGA/BKf9ti971YZ+z0z8K/+H8XF297Cp0Tr/bHziVwFOOoU4cw5gy7Vsjv+38n+/VsR3nbixtJQkPb0HAba0pzxbhoU6cN6/fanrfE+1vI+1bfx5rR/yDRlUjz+r5Fpvzn9I3C8jHjtlNJam7nJqVk+CPXdm9iBxnP2/Y1dVPiuf0SvyP9ls73k7fpIk5tEDpTfkifIq4+d0HxfkL9WYCbmDi/kkpq9g7Hdvo1pG5JjO83vnjbpzSTM/zXfEDGgOLt2Np+U9W0O1qzcJt1Jj918rv0aWhjtAX72GJqO32r1j/jZ+q6qcXWiPxtg/DkH9ovWFmmpCnAEO9Io0ysk3muMWYzsE9ETvD6D64BJgfU8Ua44mJghinjWNpvLvuGxdf4L+hxb/oX53jq9KcAeOwUa5++8tMri1+up7c4nX6t+rHwqoV8c9k3TL7Ar+s8xuMYGvdon0fp3LSO47zXf3k9o78bzRMLnE6uN7fZBTGKjNX+F7a7sDjvvI/PY+LyieQW5bJhz4bi9MkXTOaFM1/g4o6DmHNFqBM9PcEOMxsydQhdJ3YFoFlqM+K80+TfOM9vF5+23t+dfrrvIUJvlMBBj/9hzc25CoADG45MFNeCnaeUmLfv5/ATlzxFoWPRfeSH6e1c39k5HyO51Uu8+ZdjeWW5fWle1M6uD/D+cP/8iZjE33liub+1tXToUuaOPoNasdYk8veZf7cyFliz0/Xd7epmJ7VtQLNUa3Zbsj18JMw7+nfgiSv9bamhPXtzWceSI+5elXUVpzS31+nbHH+PY9Txo1h+61PUircv0O1pT/Llb18U3zMAPRv7BxDecEGQY9vL6wPeprPH6dheuWslc7f4zbHXnpxJnUR73d9cWfL6B/H1vufa6dc60rLqZrFq1yqeuqIzqW1sI8Ql8bz3l9MQERqnNCY708PcUWc4HLc+jm9iFWTdtFDf4usDX+ecTDvEdGfeTgo9hYxbPA6wynRYz8Es+sfdPD3w9pC6N3395xBrQ1rWaIdZC2DFHmvG+m7Id7j2n0TBrhM4+JvzN/ro3bQ3Zzf7syPNleT3t8S54ph0njU5JTW11oPv7urLpVP9C+v0bdOJ87NsIy+l9VOOYz3wo+3B+hS+j30rw68fs/y+EoI3BlDR4aqDRSQH6A18IiLTAYwxy4H3gBXAZ8BNxhif+r8BeAnrkF4L+N5iLwP1RGQN8A+ghNVvnHSq14n0xHRiJIZL2zuXvHyqr/8CntnS30V/e6XtKv+tx98A+8ekJ6bTunbr4jK+iTQAc6+ci8s7tX7JNUs4r815xXlT1vqDd93a41bH+S/vaHseY3qP4ZiG/jHK/5r/L3q92YtBHw8qTmtduzUnNrMtyOS4ZJ47w6/1nzvjOd4fFGpuerLvk4799861NmVfb+n0FqcTF4EQBkX77OxUT26rQ5SsOIW7exQ7F8OxcNQFfHaRfzjvTdk3Me/KeSwc+i1TLggfSK0gTG/n5p5Xhyjgexf4u/4x3tADTQNmq6Zk+pXs46dYO3/DtETuOdGvZHxKAey96aNrfftifvqnpx31fVzUK5WrP7cNiIkDJiIi1E0s2Qx3W8/baJxiBwrMzJlZnC4ixLlimHmpTdu4b2NYc4Tv3nljtb81G9i7zW7UhfeG9WFC/wlc28m/fOeo72yAvEf6WKVxcjPbY3p7ld/8NOOSGUw802mOnL/V+oh8z0TXBl0xGIfte3y/5+mVaX9zu/R2/LL7FxrWSqR2A6tMpw6eyoQ/9eTnBweS611e8+uNoQNDutbv6ghvsWCrvwcwopd9raQkxCIiLB26lNlDfuKT80KP4zP1+Qj8P33UTqjNopvG8edOt3FVt34suWYJ31z2DR+f/zHt09vz9jn2utx76p9D6gJ8fak9b+s6/nfPpQO/YfZ2v7/l0wvtoBZfQ8CH73n30bNxz+IybdK6sOy+s4t7lWD/8w1jzyEl4dDDvys6KmmSMaa5MSbBGNPIGNM/IO8hY0wbY0wHY8y0gPT5XlNUG2PMzb5egTEmzxhziTGmrTGmlzEmfFMmCAlwTdzT+x7apbfjuMbH8c1l39C3Zd/ivDhXXPFFm7zW9gxa1Qp90f37VKcJ5YqOVzgeGBHhoZMfcvRQwCqF67peR1q8v+X692P/XlzntYGv8e654f0ab539Vkhan+Z9mHvlXL4b8h19mvehUUoj7j/x/uL87y//ng51nWELsuo5V+wa1GYQ5SE3JzRE9/XdbqCsAbgWHyI+kqfo0JPU8jZfBCaeA2v/TtH+dhQd8Ecfvb7b9dRNiadZajOWDl3K0qFLGd59OImxiSTEusisncmPVzgdkv+74H9sGHtOWJNhclwy31/ud8Buyw0d3gn+BzOQAZn+7n//DH8L7M5vQmeGg/9eW7bTOpp7NbH+qp6NbOv9kv/5l28NbET8eMWPTDpvEouvWcxPV1u79o3dbyQ2JvThHtTa/3/7GgRFHmeL2qcMT2/p9E19eN6HzLliDlMHT3WYXo9rfBx/P+5mXh9og9NtO2iv0bk+M1e88z/9YNAHNEhuQI9mGfx01WLeP2eSI/+mbDuwIaNWBgB3fetfJMvnMwFoV6cdv+79lYOFB8n1rvrWqlYrTu9oQ3H4/k+fv+ac1ufw8MkP89rA14obcVdm2fv5L5/7Z9GnJ4ZGuk1NiKVlev1ipQGQHJtcbOrzMb7feEYf7x+occ8J/t7rHf07ct/5XRAR0hPTaVOnDR+e9yFd6ncBICUuhSubvEjRgbbEHMzmnyeMYfpF06mfZM18MRLDP3tbf+K0DZ9y/2z7rF/W4TJapNmeZOB7CEKfd7CN4Vf6v8LHF75NakIsY/vY4cn1EuuF/OelEXVhtz867yMm9J8Q9ga4obvTFOIKE2CrXytn3JmRx4dfiCRGYlg6dCkzL53JbcfexrVdbDfyh8t/YOrgqbzS/5WQP7JTvU7MuGRGsVkBYPL5k+naoCvhSIpNonaC3847uN1gFly1gIVXLyw2FZSG76b8ekMpS5kG0Dw+NOLsua3PLp75mZYYS8OkRmHrXr1n7yFvpoPrby41v2h/B3wRljwFjcjdeB25v11P/rYBHNx4DTcfU3p9sA/ggqsW8Mzpz7B06FIyamcAzv/108H+F32t+FrFrWsfwfu+B9NHoIIGHK3T7363c158pigfV2Vd5dj3mQafO9P2DHfn7w57vJS4FNqmtyVGYoiNiWXp0KXckO2/j30vWoC7ejlXIvQNlQbba1109SKSA8I5+5yew7oNo316e0QkbGMJ/PfSoQhsrMS6YuhYv23xC/S5M54rvp9bpjlnyx/T8BjH81joDTU95ocxYc/TMNmOFvziVzt44ORmJzOozSCHUr3rOOf1+EvXv1AaV3S8grS4NI5tdCxzrrQKNLBXWTuhNpd1uIynT3+acWeO49IOl5Z0qLCMOOsEVt44icU3vM7FHS6iaarT1n9x+4tD6ow6fpRjf84VcxjTewwLr7a+h1f6v8LxjY9n3pXzAPtOCzQV+npFMy+bWS5Za3xIDEd01TJw7/f3smHvBl7u/3KxfT6YHzf/yH0/3MfUwVPDKo/qzEtLX+LJhdZM4Gv5/WnUg1wa9xmjXCNwtwqv6AY2+hujTx1K9/umFQ9fPbPlmTzR1/pPFvz6B20apPDqyud4eVloWOnHt+1gwIGDYcNj+LBhQYqYcEMtxk2PYeH2H3AfzASE+PpfUbDjdK7u1YXaSXEM6NKYS8fN5kCBmwuym7I7t5BX/69XiccuCwu2LqB17dZhGw293uxV3DINbDH7CLTRh8vfsGeDwzT4/eXfhyjvkd+OZOq6qUy5YAqZtf09ocBjn9L8FJ49I3QeRWn88scv/Lb3N85oFbpM5Z+n/5n+mf25pP0lYWqC2+Mu8z3+5so3GTt3LJPPn+wwfczcOJNbZtzCpe0v5Z7eoT4gYwxFnqIQs2bg715yzRKHQty0fxP9P/T3xO4/8X4GtxtcYv3F1ywuNv8FMvDDgeTszwl7jrKyJ38P8a74kIZeZfHUwqd4cemLDMwYyGOnPnboCodJaSExjjrFcDSwbvc6Gqc0Lm4d+kIUHNOyDu9efxzvrvqAx+Y/gvHE8v2QOYybtYLb+nWj0G3Iuvczjsk0TBh6UtgXaKGnkB6vh4aMXrL+N4TwcZN8+OJFBceUCiQwbs7evEIO5BfRpHblP5D57nwmr5nMma3ODGvbL3QXcvus27mn9z3F3f9gDqU8SiKw3tTBU0tstUcbRZ4iVu1aFbY3UuQp4pjX/a3/74Z85+g9A/R9ry87cu3chJKutzGGtbvX0ja9bQQljw6if2lPxUHrOq0dJoOXvJOPXrjqWOJd8VzU7hLcB1uS9/vl1E6K564B2cS6YkiKd/He9b159Zr+YZUCENLL6p2byzu/by72QqS5w4fmcOc1CUk7LiP8OXzUSow7IkoBIMGVwKUdLi3R4RvniuPJ058sUSkATOhvV1ELdI6Xhe+GWPNTekL6UaMUAGJjYks0UQX7UIKVAsCXF9uh3L7rHg4RUaVwGGiP4Sjl9dkbOKV9A1rVC7PU4CFwmAC8PYUVnlZ0ivmVTbEu+rcInRG9/+e7MW7rqPT1Cg4WFPHL1v1MW7aF9Tv2c0r7Blx5/NHzYlRKZ0fuDt5e9TY3Z998WCYgpXSOnoV6lDJzddAiMuXhi22Gfg3tg+p7XN91n8Z9MRNpXOTm+Nw85iT512jwFKZh3FYB/fUM/ypSyfGxdG9Rh+4t6hy2LEr0Uj+pPrccc0tVi3FUoqYkpdwsPv4FrlrX1TGh53W3HfUTA7y0xTns88CaUfhUyC2na7deUao7qhiUcnMgtRXP51+J+0D74rRgg2TzQjvc8NtfcyjvQuSKolQtakpSIoIJevm/tWkrq+LjqONxOqNdaitWlGqP9hiUSiHd46F3mKU8g9elVRSl+qGKQSk3wWvLKooSXahiUMpNVpNa/DgydJZtSVx7UuahCymKUm1QxaAcFo1rJ7Jh7Dl0zxtPdt4LgPDngnCLisC9g0KjUiqKUn1R57NSIfbgj67p0dFHihIVaI9BiRjiHbT6tbs7nfJKDlOgKEr1RhWDEnE8xHCQxEMXVBSlWlLRFdweF5FVIrJERCaJSJ2AvJEiskZEVotI/4D0Y0VkqTfvKe8Sn3iXAX3Xmz5HRDIqIpty5Jnl6c7bRX1Z2PVe6qXEV7U4iqIcJhXtMXwBdDHGdAN+BkYCiEgnYAjQGRgAPCcivqDvzwPDsOtAt/PmA1wH/GGMaQs8ATxaQdmUI0wRsYws+gt3XHo6C+5xLnj0zBXH8MlfT64iyRRFKQ8VXdrzc2OMb+3AH4Hm3u3zgXeMMfnGmPXY9Z17iUgToJYxZrZ3Sc/XgAsC6kz0bn8AnCEaUjFqOLdbUzo3DQ2drChK9SOSPoZrAd/azs2AjQF5Od60Zt7t4HRHHa+y2QPUC3ciERkmIvNFZP727dsj9gMURVGUMgxXFZEvgcZhskYbYyZ7y4wGioA3fdXClDelpJdWJzTRmPHAeLDrMZQovKIoilJuDqkYjDFnlpYvIkOBc4EzjH/VnxwgcAX15sAmb3rzMOmBdXJEJBaoDewqw29QqpDHLu7GnR8sqWoxFEWJIBUdlTQAuAs4zxhzMCBrCjDEO9IoE+tknmuM2QzsE5ETvP6Da4DJAXWGercvBmaYmr683FHApT1bsOqBATw0OMwSjWEWZ1cUpfpT0ZnPzwAJwBdeP/GPxpjhxpjlIvIesAJrYrrJGOP21rkBeBVIwvokfH6Jl4HXRWQNtqcwpIKyKUeIxDgX6clhhqfethry9h55gRRFqRAVUgzeoaUl5T0EPBQmfT4Q0rw0xuQBl1REHqXqiHeF6R2kNrQfRVFqFNrXVxRFURyoYlAiQtuGqYcupChKjUAVgxIRMuqnVLUIiqJECFUMiqIoigNVDIqiKIoDXahHiRif//0UdOaJotR8VDEoEaN9o7SqFkFRlAigpiRFURTFgSoGRVEUxYHU9HBEIrIPWF3VchwG9YEdVS3EYVJTZVe5jyw1VW6oubKXR+5WxpgG4TKiwcew2hjTs6qFKC8iMr8myg01V3aV+8hSU+WGmit7pORWU5KiKIriQBWDoiiK4iAaFMP4qhbgMKmpckPNlV3lPrLUVLmh5soeEblrvPNZURRFiSzR0GNQFEVRIogqBkVRFMWBKgZFURTFgSoGRVEUxYEqBkVRFMWBKgZFURTFgSoGRVEUxYEqBkVRFMWBKgZFURTFgSoGRVEUxYEqBkVRFMVBxBWDiEwQkW0isqyEfBGRp0RkjYgsEZEeAXkDRGS1N29EpGVTFEVRDk1l9BheBQaUkj8QaOf9DAOeBxARF/CsN78TcLmIdKoE+RRFUZRSiLhiMMZ8A+wqpcj5wGvG8iNQR0SaAL2ANcaYdcaYAuAdb1lFURTlCFIVS3s2AzYG7Od408KlHx/uACIyDNvbICUl5diOHTtWjqSKoihRyoIFC3ZUpzWfJUyaKSU9NNGY8XgXpOjZs6eZP39+5KRTFEU5ChCRX0vKqwrFkAO0CNhvDmwC4ktIVxRFUY4gVTFcdQpwjXd00gnAHmPMZmAe0E5EMkUkHhjiLatUELfHoCv1KYpSVipjuOrbwGygg4jkiMh1IjJcRIZ7i3wKrAPWAC8CNwIYY4qAm4HpwErgPWPM8kjLd7Th9hjajPqUhz5ZWdWiKIpSQ4i4KckYc/kh8g1wUwl5n2IVhxIh3B7bU5jw/XruPldH/ypKJCksLCQnJ4e8vLyqFqVEEhMTad68OXFxcWWuUxU+BuUIYrz+e49akhQl4uTk5JCWlkZGRgYi4cbPVC3GGHbu3ElOTg6ZmZllrqchMaIcdS0oSuWRl5dHvXr1qqVSABAR6tWrV+4ejSqGKMejmkFRKpXqqhR8HI58qhiiHDUhKYpSXlQxRDnaY1CU6MblcpGdnU2XLl245JJLOHjwYIWPqYohylG9oCjRTVJSEosWLWLZsmXEx8czbty4Ch9TRyVFOTqxTVGODPf9bzkrNu2N6DE7Na3FPwd1LnP5Pn36sGTJkgqfV3sMUY76GBTl6KCoqIhp06bRtWvXCh9LewxRjvoYFOXIUJ6WfSTJzc0lOzsbsD2G6667rsLHVMUQ5ahiUJToxudjiCRqSopyVC8oilJeVDFEOdpjUBSlvKhiiHLU+awo0c3+/fsjfkxVDFGODldVFKW8qGKIclQvKIpSXlQxRDnqY1AUpbxUxgpuA0RktYisEZERYfLvEJFF3s8yEXGLSF1v3gYRWerNmx9p2Y5G1MegKEp5ieg8BhFxAc8C/YAcYJ6ITDHGrPCVMcY8DjzuLT8I+LsxZlfAYfoaY3ZEUq6jGe0xKIpSXiLdY+gFrDHGrDPGFADvAOeXUv5y4O0Iy6AEoM5nRVHKS6QVQzNgY8B+jjctBBFJBgYAHwYkG+BzEVkgIsNKOomIDBOR+SIyf/v27REQO3pRU5KiRDcPPfQQnTt3plu3bmRnZzNnzpwKHzPSITHCLRVU0qtpEPB9kBnpJGPMJhFpCHwhIquMMd+EHNCY8cB4gJ49e+qrrxS0w6Ao0cvs2bOZOnUqCxcuJCEhgR07dlBQUFDh40ZaMeQALQL2mwObSig7hCAzkjFmk/d7m4hMwpqmQhSDUnbUx6AoR4hpI2DL0sges3FXGDi2xOzNmzdTv359EhISAKhfv35EThtpU9I8oJ2IZIpIPPblPyW4kIjUBk4FJgekpYhImm8bOAtYFmH5jjpUMShK9HLWWWexceNG2rdvz4033sisWbMictyI9hiMMUUicjMwHXABE4wxy0VkuDfft7TQYOBzY8yBgOqNgEnehatjgbeMMZ9FUr6jEdULinKEKKVlX1mkpqayYMECvv32W77++msuu+wyxo4dy5/+9KcKHTfiYbeNMZ8CnwaljQvafxV4NShtHdA90vIc7WiPQVGiG5fLxWmnncZpp51G165dmThxYoUVg858jnJ0VJKiRC+rV6/ml19+Kd5ftGgRrVq1qvBxdaGeKEd7DIoSvezfv59bbrmF3bt3ExsbS9u2bRk/fnyFj6uKIcrRCW6KEr0ce+yx/PDDDxE/rpqSohzVC4qilBdVDFFOxH0MBQfhyzFQmBvhAyuKUl1QxRDlRNzHMP9l+O4J+OHpyB5XUWoo1d1cezjyqWKIciKuGHzHy90d2eMqSg0kMTGRnTt3VlvlYIxh586dJCYmlqueOp+jnIjfr64474HdET6wotQ8mjdvTk5ODtU5mGdiYiLNmzcvVx1VDFFOxHsMMd5bxlMU2eMqSg0kLi6OzMzMqhYj4qgpKcqJuPM5xuU9sCoGRYlWVDFEORHvMYgqBkWJdlQxRDsR7zH4TEnqY1CUaEUVQ5RTaSExVDEoStSiiiHKibiPodiEVD2H5ymKUnFUMUQ5Ee8x+BRDNR23rShKxYm4YhCRASKyWkTWiMiIMPmnicgeEVnk/dxb1rpK+Yn4xBufCcl4IntcRVGqDRGdxyAiLuBZoB92/ed5IjLFGLMiqOi3xphzD7OuUg4qzZSkikFRopZI9xh6AWuMMeuMMQXAO8D5R6CuUgKVZkpSH4OiRC2RVgzNgI0B+znetGB6i8hiEZkmIp3LWRcRGSYi80VkfnWeil4d0B6DoijlJdKKQcKkBb+aFgKtjDHdgaeBj8tR1yYaM94Y09MY07NBgwaHK+tRQeX5GLTHoCjRSqQVQw7QImC/ObApsIAxZq8xZr93+1MgTkTql6WuUn4i/v7WUUmKEvVEWjHMA9qJSKaIxANDgCmBBUSksYiId7uXV4adZamrlJ9K8zF4CiN7XEVRqg0RHZVkjCkSkZuB6YALmGCMWS4iw73544CLgRtEpAjIBYYYa+8IWzeS8h2NRNzH4Au37S6I8IEVRakuRDzsttc89GlQ2riA7WeAZ8paV6kYke8x+BSDBtFTlGhFZz5HORF3PvtGI2mPQVGiFlUMUU7kh6uqKUlRoh1VDFFO5Nd89ioGXY9BUaIWVQxRTuSHq2qPQVGiHVUMUU7kfQw+xaDDVRUlWlHFEOVUno9BFYOiRCuqGKKcyPsYvMfTCW6KErWoYohydIKboijlRRVDlFNpQfR0gpuiRC2qGKKcShuuqj0GRYlaVDFEOTpcVVGU8qKKIcqpNB8Dxq8kFEWJKlQxRDmRD6IXsHKb9hoUJSpRxRDlVNoEN9C5DIoSpahiiHIib0oK7DGoYlCUaEQVQ5QTaEqKSO8h0K+gk9wUJSqJuGIQkQEislpE1ojIiDD5V4rIEu/nBxHpHpC3QUSWisgiEZkfadmORjwBXYaI9B4cpiT1MShKNBLRFdxExAU8C/QDcoB5IjLFGLMioNh64FRjzB8iMhAYDxwfkN/XGLMjknIdzbhNoGIwuJCKHdCjPgZFiXYi3WPoBawxxqwzxhQA7wDnBxYwxvxgjPnDu/sj0DzCMigBuANcAhHxQ6uPQVGinkgrhmbAxoD9HG9aSVwHTAvYN8DnIrJARIaVVElEhonIfBGZv3379goJHO14gnoMFT+g+hgUJdqJqCkJwtopwr6NRKQvVjGcHJB8kjFmk4g0BL4QkVXGmG9CDmjMeKwJip49e0Z63E1U4fYEOp8jcED1MShK1BPpHkMO0CJgvzmwKbiQiHQDXgLON8bs9KUbYzZ5v7cBk7CmKaUCVEqPQVx2W01JihKVRFoxzAPaiUimiMQDQ4ApgQVEpCXwEXC1MebngPQUEUnzbQNnAcsiLN9Rh3NUUgQUg/FAbKLdVsWgKFFJRE1JxpgiEbkZmA64gAnGmOUiMtybPw64F6gHPCciAEXGmJ5AI2CSNy0WeMsY81kk5TsacY5KisABjRtiE6DwgJqSFCVKibSPAWPMp8CnQWnjArb/DPw5TL11QPfgdKViBI5KCu/tKSceD8QlQS7g0TUZFCUa0ZnPUU7kTUneHgNoj0FRohRVDFFO8AS3CuNxq49BUaIcVQxRjifSPgZPEcQl221VDIoSlahiiHI8nggH0XMXQnyK9+CqGBQlGlHFEOW4A3RBRGYCugv8ikF9DIoSlahiiHIi7nz2FKopSVGiHFUMUY470mG33YUQr4pBUaIZVQxRjmNUUiQ0g7sA4tSUpCjRjCqGKMcTySB6xjh7DDrBTVGiElUMUU5Eg+h53ICx8xgkRnsMihKlqGKIciI6Ksk3PNUVB6549TEoSpSiiiHKCTQluSvqY/D1EGLi7EcVg6JEJaoYopxAZVDoiKh3OAfz9Rjiba9BJ7gpSlSiiiHKCRyVlF9UUcXg7TG4Yq1iUB+DokQlqhiinIIiD3Euu+JqfqH7EKUPdbAD9js+TX0MihLFqGKIcvKLPNRKjCverhB5e+13QhrExKpiUKonezdBkfZmK0LEFYOIDBCR1SKyRkRGhMkXEXnKm79ERHqUta5SfvKL3NRKipBiyN9jvxNreXsM+vAph89t7y0mY8QnXPz8DxU7kDF2KPVvc2DxO/CfLHgqOyIyHq1EdAU3EXEBzwL9gBxgnohMMcasCCg2EGjn/RwPPA8cX8a6FccY78cDGGsrrwJ8kU59LgCfJ8B4PAH7xuYbg8E49gmbH3hcgzGG/fv20q5RKlt25FGUtw/yU3Ce0Xm8Uve3LrffqY0gLhF2roHlk6DVyZDaoKQfaq+1x22/jQfjKbLf7iJMwX6ITcK4CzAHd2FikyApHZOQijFiHdzuQkxRASAggrhiMQgeBCMuXDbZ5uP7b50IHjAQI0CMC4/Hg8dt8BgPHmMQiSHWJYAQExODiCAiYAxi3N4TCCIxdg6H2DaVTfanSUwJba3g+8677fF4cBsPbrfbymQM7iIP4CFGwCUGMQT8Lm9dDEhssRwYt/8aeGU1CAaDxwgeDxhvuk9wiRFcYsubgHRj8KfZg9r/QgSKvw3GUwieIozHjXiK7P9k7LZxF4LHTYxxI55CYnBTVFRInHjIy89n44597Fy0gmtc29i1sRYLVzaiZd0E3EUFuDxFJLk8GHchpqjQfrsLMAX7yc3LJ+/AXhqkp5OU+zuydgYxv0zHJNdDDu70X++9v8OG7yB/P2xaiMnfBx0GQvPjMLiKf1fx9fD+lyJCUVERhQX5iBgMMZgYe7/5njNjfE+b9xoVPyYe+7d6E1wxQlJsDG6P8d5jYBBcMUKM+L8lxhXwvxEgW9UhEQnF7DuYSG9gjDGmv3d/JIAx5pGAMi8AM40xb3v3VwOnARmHqhuOHk3jzHd/SUcC/ir7jffVATG+bwn9rbkm3v9QBBFc2n8rhZ7Lt2/CpJdWJ5xM1ZnNph5neJ7hWpnC7TFvA+AxQiGxuIkhBo/3Y+ynhv2+iuIxvruu+LVx1F2DqmCpJ4OmspOHi65kqWnNh3H/JE1yq1qsCuMxgpsYr5J3vqds86BkJSIhbzDLuoQsOo36FhFZYIzpGa5MpNd8bgZsDNjPwfYKDlWmWRnrAiAiw4Bh3t38lPu3L6uAzFVFfWBHVQtRfvYC59S/A3bcUdWilJ8aes1V7kOzxPv9GAC1Kn7AKL7m38FoAWhVUolIK4Zw6qukhndwmbLUtYnGjAfGA4jI/JK0XnWmpsoNNVd2lfvIUlPlhpore6TkjrRiyAFaBOw3BzaVsUx8GeoqiqIolUykRyXNA9qJSKaIxANDgClBZaYA13hHJ50A7DHGbC5jXUVRFKWSiWiPwRhTJCI3A9MBFzDBGLNcRIZ788cBnwJnA2uAg8D/lVa3DKcdH8nfcASpqXJDzZVd5T6y1FS5oebKHhG5IzoqSVEURan56MxnRVEUxYEqBkVRFMVBjVYMNSWEhoi0EJGvRWSliCwXkb9508eIyO8issj7ObuqZQ1GRDaIyFKvfPO9aXVF5AsR+cX7nV7VcgYiIh0CrukiEdkrIrdW1+stIhNEZJuILAtIK/Eai8hI7z2/WkT6V43UJcr9uIis8oa7mSQidbzpGSKSG3Dtx1UzuUu8N6r59X43QOYNIrLIm16x6+0LnVDTPlgH9VqgNXao62KgU1XLVYKsTYAe3u004GegEzAGuL2q5TuE7BuA+kFpjwEjvNsjgEerWs5D3CdbsJN5quX1Bk4BegDLDnWNvffNYiAByPQ+A65qJPdZQKx3+9EAuTMCy1XD6x323qju1zso/9/AvZG43jW5x9ALWGOMWWeMKQDeAc6vYpnCYozZbIxZ6N3eB6zEzvSuqZwPTPRuTwQuqDpRDskZwFpjzK9VLUhJGGO+AXYFJZd0jc8H3jHG5Btj1mNH9/U6EnIGE05uY8znxpgi7+6P2PlI1YoSrndJVOvr7UNEBLgUeDsS56rJiqGk0BrVGhHJAI4B5niTbvZ2uydUN5OMFwN8LiILvKFIABoZO/cE73fDKpPu0AzB+bBU9+vto6RrXJPu+2uBaQH7mSLyk4jMEpE+VSVUKYS7N2rK9e4DbDXG/BKQdtjXuyYrhjKH0KguiEgq8CFwqzFmLzaybBsgG9iM7QpWN04yxvTARsW9SUROqWqByop3ouR5wPvepJpwvQ9FjbjvRWQ0UAS86U3aDLQ0xhwD/AN4S0QiENIoYpR0b9SI6w1cjrMBVKHrXZMVQ1nCb1QbRCQOqxTeNMZ8BGCM2WqMcRtjPMCLVFEXtTSMMZu839uASVgZt4pIEwDv97aqk7BUBgILjTFboWZc7wBKusbV/r4XkaHAucCVxmvw9ppidnq3F2Bt9e2rTkonpdwbNeF6xwIXAu/60ip6vWuyYqgxITS89r+XgZXGmP8EpDcJKDYYqFZRYkUkRUTSfNtYx+Iy7HUe6i02FJhcNRIeEkcrqrpf7yBKusZTgCEikiAimdh1TeZWgXxhEZEBwF3AecaYgwHpDcSuuYKItMbKva5qpAyllHujWl9vL2cCq4wxOb6ECl/vqvCuR9BLfzZ2hM9aYHRVy1OKnCdju59LgEXez9nA68BSb/oUoElVyxokd2vsiIzFwHLfNQbqAV8Bv3i/61a1rGFkTwZ2ArUD0qrl9cYqr81AIbaFel1p1xgY7b3nVwMDq5nca7A2ed99Ps5b9iLvPbQYWAgMqmZyl3hvVOfr7U1/FRgeVLZC11tDYiiKoigOarIpSVEURakEVDEoiqIoDlQxKIqiKA5UMSiKoigOVDEoiqIoDlQxKIoXEakXEI1yS0C0zf0i8lwlnfNWEbmmlPxzReS+yji3opSEDldVlDCIyBhgvzHmX5V4jljsGPMexh94LriMeMucZAImjClKZaI9BkU5BCJymohM9W6PEZGJIvK5N/79hSLymNg1Kz7zhj5BRI71Bi9bICLTg2bW+jgdG7KjyFvnryKywhvI7R0AY1tuM7EhJhTliKCKQVHKTxvgHGxI5jeAr40xXYFc4ByvcngauNgYcywwAXgozHFOAhYE7I8AjjHGdAOGB6TPx0bPVJQjQmxVC6AoNZBpxphCEVmKXQjoM2/6UuwCKR2ALsAX1hKECxvKIJgm2LU5fCwB3hSRj4GPA9K3AU0jJ76ilI4qBkUpP/kAxhiPiBQav6POg32mBFhujOl9iOPkAokB++dgV+k6D7hHRDp7zUyJ3rKKckRQU5KiRJ7VQAMR6Q025LqIdA5TbiXQ1lsmBmhhjPkauBOoA6R6y7WnekeCVaIMVQyKEmGMXWr2YuBREVmMjTJ6Ypii07A9BLDmpje85qmfgCeMMbu9eX2BTypTZkUJRIerKkoVIiKTgDuNc0nGwPxGwFvGmDOOrGTK0YwqBkWpQkSkA3Z9529KyD8OKDTGLDqigilHNaoYFEVRFAfqY1AURVEcqGJQFEVRHKhiUBRFURyoYlAURVEcqGJQFEVRHPw/j+dLYlIMJhEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "response = picker.annotate(stream)\n",
    "# Plot the response\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "time_diff = response[0].stats.starttime-stream[0].stats.starttime\n",
    "fig, ax = plt.subplots(2, 1, sharex=True)\n",
    "t0 = np.arange(stream[0].stats.npts)/sample_rate\n",
    "t1 = np.arange(response[0].stats.npts)/sample_rate + time_diff\n",
    "ax[0].plot(t0, np.array([x.data for x in stream[:3]]).transpose(), label=[x.stats.channel for x in stream[:3]])\n",
    "ax[1].plot(t1, np.array([x.data for x in response[1:3]]).transpose(), label=['P', 'S'])\n",
    "ax[1].set_xlabel('Time (s)')\n",
    "ax[1].set_ylim(0, 1)\n",
    "for j in ax:\n",
    "    j.legend()\n",
    "    j.set_xlim(0, t0.max()+1)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3077d818",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Retrieve the predicted picks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4e370f53",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Printing picks\n",
      "CN.BJ.\t2018-07-18T15:25:01.800000Z\tS\n",
      "CN.SC.\t2010-03-14T20:35:00.900000Z\tS\n",
      "CN.BJ.\t2018-07-18T15:24:51.120000Z\tP\n",
      "CN.SC.\t2010-03-14T20:34:55.260000Z\tP\n",
      "{'P': {'trace': ['CN.BJ.', 'CN.SC.'], 'time': ['20180718T15:24:51.12', '20100314T20:34:55.26']}, 'S': {'trace': ['CN.BJ.', 'CN.SC.'], 'time': ['20180718T15:25:01.80', '20100314T20:35:00.90']}}\n"
     ]
    }
   ],
   "source": [
    "classifyoutput = picker.classify(stream, P_threshold=.3, S_threshold=.3)\n",
    "picks = classifyoutput.picks\n",
    "print('Printing picks')\n",
    "lines = {phase:{'trace':[],'time':[]} for phase in 'PS'}\n",
    "def format_utc(t):\n",
    "    return '%04d%02d%02dT%02d:%02d:%02d.%02d'%(t.year,t.month,t.day,t.hour,t.minute,t.second,t.microsecond/1e4)\n",
    "for pick in picks:\n",
    "    print(pick)\n",
    "    # Please ensure both the station and network code exist in the header of the SAC files\n",
    "    lines[pick.phase]['trace'].append(pick.trace_id)\n",
    "    lines[pick.phase]['time'].append(format_utc(pick.peak_time))\n",
    "print(lines)\n",
    "import json\n",
    "# save results as json\n",
    "with open('results/picks_demo.json','w') as fp:\n",
    "    json.dump(lines,fp)\n",
    "# read results as dict\n",
    "with open('results/picks_demo.json','r') as fp:\n",
    "    data = json.load(fp)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
